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ABSTRACT 


A Galerkin-Weighted Residuals fomulatlon is employed to establish 
an implicit finite element solution algorithm for generally non-linear 
initial -boundary value problems. Solution accuracy, and convergence rate 
with discretization refinement, are quantized in several error norms, by 
a systematic study of numerical solutions to several non-linear parabolic 
and a hyperbolic partial differential equation characteristic of the 
equations governing fluid flows. Solutions are generated using selective 
linear, quadratic and cubic basis functions. Richardson extrapolation is 
employed to generate a higher-order accurate solution to facilitate 
isolation of truncation error in all norms. Extension of the mathematical 
theory underlying accuracy and convergence concepts for linear elliptic 
equations is predicted for equations characteristic of laminar and turbulent 
fluid flows at non-modest Reynolds number- The non-diagonal initial -value 
matrix structure introduced by the finite element theory is determined 
intrinsic to improved solution accuracy and convergence. As an alternative 
to the conventional multi -dimensional finite element algorithm, a factored 
Jacobian iteration algorithm is derived and evaluated to yield a conse- 
quential reduction in both computer storage and execution CPU requirements 
while retaining solution accuracy. The developed hypermatrix statement of 
the solution algorithm reduces storage requirements and facilitates direct 
inclusion of parameter variations. The results of the research conducted 
under the Grant and reported herein document an accurate and versatile 
algorithm potentially applicable to solution of a wide range of practical 
problem classes in aerodynamics and fluid mechanics. 


INTRODUCTION AND SUMMARY 


The primary objective of this research project i*s to assess predic- 
tion of extension of the mathematical theory governing accuracy and 
convergence character of finite element solution of linear elliptic partial 
differential equations, to the progressively more complex hyperbolic and 
non-linear parabolic partial differential equations characteristic of fluid 
mechanics. The results for the non-linear laminar and turbulent flow 
cases, cofnsidered and reported herein, predict that extension of the linear 
equation theory Is valid for the finite element solution algorithm using 
linear, quadratic and cubic elements. Comparison tests between the finite 
element and a finite difference (Crank-Nicholson) algorithm have quantized 
for the first time the differences in numerical accuracy attainable. These 
comparison results generally confirm that the linear finite element solution 
algorithm is consistently superior to the equal order accurate Crank- 
Nicholson algorithm in terms of accuracy and convergence, while maintaining 
comparable solution economy, for the non-linear parabolic equations 
considered. The reported results of the tightly controlled numerical 
experiments confirm viability of the energy norm as the intrinsic measure 
for accuracy and convergence determination in laminar and turbulent 
parabolic-type flowfield prediction. This is in significant distinction to 
the variability in convergence measured in the various common engineering 
norms. 

A range of practically useful finite element discretizations for 
parabolic flow prediction has been observed. Solutions employing coarse 
grid linear finite element discretizations generally display accuracy 
superior to those predicted by strict adherence to the convergence curve. 
Conversely, those solutions obtained with quadratic finite elements typically 
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display coarse grid inaccuracy. On the other hand, on progressively refined 
discretizations, sources of error other than that associated v.'ith finite 
element discretization serve to obliterate the refined solution accuracy 
theoretically obtainable. The absolute error associated with solutions 
obtained employing quadratic elements is, however, uniformly smaller than 
that associated with linear element solutions on sufficiently refined grids. 
The accuracy obtainable using a non-uniform grid within the finite element 
algorithm was found to be superior to use of uniform grids for the parabolic 
problems studied. This is not uniformly true, however, for the finite 
difference algorithm evaluated. Furthermore, while non-uniform discretiza- 
tions display better absolute error for both the linear and quadratic finite 
element algorithms, the presence of an optimum-accuracy grid was detected 
for linear element solutions, but was absent for quadratic element solutions. 

The use of non-zero pressure gradients for the laminar parabolic 
flows did not measurably alter the level of accuracy or convergence character, 
measured in the energy norm, of the linear finite element algorithm. This 
was not the case, however, for the quadratic finite element algorithm, where 
the fourth order accuracy of the algorithm for zero pressure gradient was 
degraded to second order for the cases involving non-zero pressure gradient. 
This may be due in part to the alteration in the convergence character, from 
oscillatory for zero pressure gradient, to monotonic for the non-zero 
pressure gradient cases. Computation of the transverse velocity distribution, 
using generally second order finite difference formulae, yields a significant 
source of error in actual computations on fine grids, which adversely affects 
the accuracy attainable using higher order accurate finite element interpola- 
tions. The level of error, capable of quantization, increased from 10"® to 
about 10”®, when the non-linearly induced error, stemming from the transverse 
velocity solution methodology, was removed. Hence, while transverse velocity 
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constitutes data within the theoretical framework of the algorithm, a uni- 
formly fourth-order accurate algorithm would be required to increase solution 
accuracy beyond about 10“^. The results of the turbulent boundary layer 
solutions indicate that a strict! y-accurate evaluation of the Jacobian, within 
the Newton iteration algorithm, is not necessary to achieve an adequately- 
accurate engineering solution. Significant solution economics can result, 
therefore, in terms of computer core and CPU, by taking advantage of the 
versatility embedded within the developed modified Newton iteration algorithm 
for multiple dependent variable systems. 

The transient continuity equation solutions confirmed that, on all 
comparison bases, the performance of the implicit finite element algorithmic 
solution form for a dominantly hyperbolic equation, is superior to the 
equivaleht-cbmplextty finite difference form with no additional computational 
effort. The primary objective with the conducted numerical experiments was 
to evaluate economy measures applicable to the basic finite element formula- 
tion, and to assess their influence on determined accuracy and convergence. 
The developed factored Jacobian integration algorithm displayed considerable 
economy in terms of computer storage and CPU, in comparison to the conven- 
tional multi -dimensional finite element algorithm, with no measurable loss 
in accuracy. Accuracy and convergence properties of the factored algorithm 
have been quantized, and the several numerical solutions obtained for a 
variety of velocity fields, document accuracy and computational aspects and 
illustrate its versatilTty. 

THEORETICAL ANALYSIS 

Accuracy and Convergence 

In finite element analysis, error estimates and convergence properties 
are typically expressed in an energy norm, cf. Strang and Fix (1973). 
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Alternatively, for finite differences, a stability analysis is employed to 
ascertain that the method is convergent, and a local truncation error analysis 
determines order-of-accuracy by means of a Taylor series expansion. Error 
and convergence may also be measured in other norms including the familiar 
engineering parameters. For the case of boundary layer solutions, for 
example, these could include the integral parameters of boundary layer 
displacement and momentum thickness, shape factor, and skin friction coeffi- 
cient. 

The primary focus of this reported analysis is numerical determination 
of the accuracy and convergence character of a finite element numerical 
solution algorithm for representative initial-value problems associated with 
high Reynolds number laminar, turbulent, and inviscid flows. A companion 
focus is evaluation of error measurement norms that facilitate estimation of 
contributing fac+nrs to solution inaccuracy, hence solution economy. The 
question of accuracy of a finite element (or any other numerical) solution 
requires quantization with reference to acceptable (usable and efficient) 
interpolation functions and mesh size distributions. The mathematical theory 
of finite elements, which examines these details in thoroughness, is 
generally limited to linear partial differential equations. The present 
requirement is to recall the fundamental theoretical concepts as applied to 
elliptic equations, and to present the extension to accuracy and convergence 
measure for the hyperbolic and non-linear parabolic equations of interest. 

The point of departure, cf. Strang and Fix (1973), is the linear 
differential equation on one dimensional space 

Uq) = V ^ + g = 0 (1) 


where v is a general (constant) diffusion coefficient, and g is a source 


( 2 ) 


term. The boundary conditions are typically assumed homogeneous as 

q{Xjt) = 0 

To establish the error measure for the finite element solution of equation 
(1) for m= l, which corresponds to elTiptic, the fundamental requirement is 
determination of how clpse the finite element solution q* 

^ "."t . , ' ^ ■ ' 

q* 5 r<is - {3) 

e'i"! ® 6^1 * ® ' 

is to the true solution q(x). The fundamental theorem states that the finite 
element solution lies as close as possible to the exact solution, in the 
sense that the energy in the error, e 

E = q - q* (4) 

is minimized: The minimum energy functional that is the equivalent to 
minimization of the variational statement for equations (1)«'{2) is 



The energy inner product for equation (1) is defined as 

ECq^q) - ||v (6) 

The proof Of the fUrtdamehtal theorem, cf. Strang end 
pp* 39-41, assures that the energy in the error is minimized finite 

element solution of equation (1) obtained using equation (3). Hence, the 
error in the energfy inner product, equation (6), for the finite elemeht 
solution, equation (3), satisfies the inequality 

ICe,s) 5 II q-q^ 1! < C |[ q }| ^ (7) 
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where is a constant independent of the measure of the largest finite 
element on R^, the exponent of is a function of the largest complete 
polynomial degree k in equation (3), 2m = 2 is the order of the elliptic 
operator, equation (1), and a=k + l is related to the required smoothness of 
the solution. 

Equation (7) states that the error in the energy inner product of the 
finite element algorithm, using a linear interpolation polynomial for example, 
goes to zero as the order A^ {or more conventionally in finite difference 
terminology, h^, i.e. , the method is second-order accurate). The error bound 
in equation (7) can be refined, for m^l and k = l for example, as 

E{e.e) iC^Alllq'll^ 

iCjA'IlFlP (8) 

and Cg are constants, A^ is the measure of the largest finite 
and |IF(| is the Lg norm of the data of the problem specification, 

IlFll - [Jg^dx]*® (9) 

Equation (9) states that the error in the energy of the linear finite element 
approximation to the true solution q decreases in proportion to the square 
of the measure of the largest finite el enfant, and that the error depends 
continuously on the data of the problem specification. Hence, as a consequence 
of the fundamental theorem, the finite element solution converges in the 
energy norm as 

E(s,s) — ► 0 as Ag — 


where Cg 
element, 
i.e.. 


0 


(10) 


l*or the class of time-depettdeht an4/or initia 
interests extetiSibh of equation (1) to a linear, eiTiptic, initial -vaTu^^ 
equation statenient yietd^^^^ 



with boundary condItiohS 

. q(S,t) = 0 • ^ (12) 

and an initial Gondftion 

q;(XvO) = q^^^^ (13) 

Oden and Reddy (1976) devetop* error estimates for fully discrete Gaierkin 
approximation of an equation of the form {11) * the components of the 
approximate spluti on scheme are defined as: q(x,t)y thfe^exact So^T^ 
equation (11); Q(x,t), the solution to the semidiscrete;^.(^^^^ 
dependence on time, t, is still assumed) Gaierkin approximation; and q*(x), 
the solution to the fi nite-differeiice Gai erkin approximati on - at 'time^ 

The following definitions for error are then introduced: 

eCXsrtAt) - q(x*nAt) - q*(x) = approximation error 

o(x,nAt) = q(XinAt) - Q(x,nAt) - semidiscrete approximation error 

T(x,nAt) = Q(x,nAt) - q*(x) =* temporal approximation error (14) 

For any choice of norm 

I!e(t)if “ Ijo + t|} < jjo|{ + 11t|| (15) 

For use of a forward difference integration algorithm, Oden and Reddy (1976) 
prove that if the components of the error are given by equation ( 14 ), then 
for a linear, el 1 i ptic , i ni tial -val Ue equati on , the error satisfies the 
inequality ■ ‘ ^ 
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|ls(nAt)l| < CjA^^^'^^""’^llq{nAt)i| + CgAtllQ^jljj^ (16) 

Here, the first term on the right hand side of equation (16) is the right 
hand side in the inequality (7). Furthemore, Cg is a constant. At is the 
time step, and 1!Q°H^ is the norm of the initial data. 

For the second case of interest, consider a first order hyperbolic 
problem statement of the form 

with initial condition 

q(0) = q^ (18) 

Here again, for use of a forward difference integration algorithm, Oden and 
Reddy (1976) prove that if the components of the error are given by equation 
(14), then for a linear hyperbolic equation, the error satisfies the 
inequality 

nAt 

]|s{nAt)|[ <. CjAg ^llq||j^^2 "** ^ 

where Cg is another constant. Unfortunately, no similar analyses for an 
initial -value, non-linear problem exists, which prompts the numerical 
experiment approach taken herein to study the convergence character of 
finite element solution of non-linear initial-value problems. 

Error Analysis 

The main emphasis in this analysis is assessment of the discretization 
error associated with use of the finite element solution procedure for non- 
linear problems. For such problems, the approach taken is evaluation of 
absolute accuracy as determined by computed solution comparisons, for 
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progressively refined di.scretT 2 ations, in the energy norm, to establish the 
convergence rate exponent on Ag. Since the fluid mechanics equations to be 
solved correspond to statements of conservation, correspondingly defined 
norms are also useful, in particular 


Pl = I f 'la'lf 

(20) 

«R" 

(21) 

The familiar engineering parameters useful for quantizing flow phenomena can 
typically be constructed from the energy and p-norms. For example, in 
boundary layer flow, shape factor and skin friction are parameters of great 
engineering significance in assessing solution acceptability, i.e., accuracy 
Shape factor H is defined as 

. H s a*0“^ 

(22) 

where S* is the boundary layer displacement thickness and 6 
thickness defined as 

is the momentum 

6* 5 |[1 - uu-^ldXg 
6 

(23) 

9 s juj^u'^tl - Uj^uj^ldXg 
■ 5 

(2+) 

where uj is the local inviscid freestream velocity, and u is the (time- 
averaged) boundary layer velocity distribution. Assuming u noh-dimensional- 
ized by uj, and for S spanning UR^ = R^, using equations (20) -(21) in the 


discrete approximation to equations {23)-(24) yields 
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6* = 6 - Pj 

(25) 

8 = pj - 2pg 

(26) 


Skin friction is an engineering measure of drag, hence the viscosity induced 
shear stress, and is defined as 



Equation (27) is also evaluable using the defined p^. 


DISCUSSION AND RESULTS 

Presented herein are the results of the numerical evaluation of 
solutions of the selected non-linear parabolic and hyperbolic partial differ- 
ential equations. Primary emphasis is on quantization of solution accuracy 
and convergence. with discretization refinement. Test cases used are the 
steady two-dimensional incompressible laminar and turbulent boundary layer 
flow and laminar and turbulent parabolic flow in a duct. For turbulent flows 
a comparison between the mixing length closure and the turbulent kinetic 
energy model has been established. Additional results are presented of the 
numerical solution of the transient continuity equation, with primary 
emphasis on solution economy, accuracy and convergence of the developed split 
Jacobian finite element algorithm. 

Parabolic Equation Solutions 

Problem Statement 

It is required to establish the two-dimensional velocity and pressure 
distributions, tt(x,y) and p(x,y), where 
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u{x,y) s u(x,y)i + v(x,y)j (28) 

Using the boundary layer order of magnitude analysis for large Reynolds 
number, cf. Schlichting (1968), the parent time-averaged steady flow Navier- 
Stokes equations in non-dimensional form are 


Lfp) = |H.+ |I= 0 
3x 3y 


( 29 ) 


US) = 


~ 3u . - 8u 3 I 

11 -IZ- + V 

ax ay ay' 


V - u"v' 
ay ^ 


+ -^= 0 

pdx 


(30) 


At the edge of the boundary layer, the viscous terms are zero, and the 
x-momentum equation (30) with the continuity equation (29) reduces to 


U 


aUj 

IIS' 


- IM 
" p ax 


(31) 


where Uj is the inviscid flow velocity at the edge of the boundary layer. 

The x-axis is assumed aligned with the direction of predominant flow, and 
y is the coordinate traversing the thickness of the boundary layer, see 
Figure 1. Under the large Reynolds number assumption, the transverse 
momentum equation is identically satisfied by a pressure distribution 
impressed uniformly across the boundary layer thickness, i.e., p(x,y) = Pj(x). 

Closure of equations (2)-(4) requires a relationship be established 
for laminar viscosity v and the dominant Reynolds stress shear component 
y'v''. Kinematic viscosity v is a property of the fluid constituting the 
boundary layer, and a constant for isoenergetic subsonic flows. Closure for 
the Reynolds stress u‘'v" is accomplished in the elementary form by assumption 
of an "effective" viscosity coefficient v® defined as 



( 32 ) 


4 * 

where v is the “turbulent kinematic viscosity" correlation coefficient 
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— s— s- _ t 9u 
■U V s V ^ 

9y 


(33) 


Hence, the parabolic equation (30) of primary interest in this investigation 
takes the form 


L(u) = u 


3x 


. ~3u 
^ dy 




(34) 


Two closure models are used in this analysis. Prandtl mixing length 
theory (MIT) establishes an algebraic relation far (cf. Schlichting, 

1968; Cebeci and Smith, 1974). Using a dimensional analysis, involves 
the product of a scale velocity and a scale length. Using the mean velocity 
gradient for the former, and Prandtl's mixing length I for the latter, 
yiel ds 


2 


30 

3y 


(35) 


The mixing length £ is defined as 



0 < y £ X5/k 
y > XS/k 


(36) 


Hie Van Driest function w accounts for the damping influence of a wall on 
the velocity fluctuations t". Following a rigorous analysis (cf. Cebeci 
and Smith, 1974), the damping function form is 

to’ s 1 - exp(-y/A) (37) 

In equation (36), y is the coordinate normal to the surface, 6 is the 
boundary layer thickness, and X and k are constants (typically 0.09 and 
0.435 respectively). In equation (37), A is a complicated function of many 
factors influencing flow phenomena near the surface including axial pressure 
gradient and normal mass flow. The form of Cebeci and Smith (1974) serves 
to unify the many formulations as 
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A s A'^vN-i 


A 'y 

-1/2 

► ^ 
0 

KJ 


• 


(38) 


where 


5:^ 


V 

Pi 


w 



+ exp 




(39) 


All variables are time-aver^ iid steady components, subscripts I and w refer 
to inviscid freestreaun and wall values respectively, is a constant (25.3), 
p"** and V'*' are functions accounting for axial pressure gradient and mass 
addition respectively, cf. Cebeci and Smith (1974). In equation (38), is 
the wall shear stress defined as 


■w 


~ Pw^w3y 


w 


(40) 


In this analysis the wall shear stress is evaluated from the Ludwieg-Tillman 
equation for skin friction 



s 0.246 X 10 


- 0 . 678 Hj^g 


-0.268 

e 


(41) 


In equation (31), H is the shape factor and Re0 is the Reynold's number based 
on the momentum thickness e. Since the freestream outside the boundary layer 
is assumed free of turbulence, the intermittency factor y in equation (35) is 


Y = [l + (y/6)®]'“^ y>S (42) 

which serves to provide a rapid decrease of at the freestream edge of 6(x). 

An alternative formulation to Prandtl's mixing length model that yields 
a differential equation statement is the turbulent kinetic energy (TKE) 
two-equation model. For this, the scale velocity is selected as the kinetic 
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energy k of the velocity fluctuations 

k s (43) 

the length scale is defined as the scale length of isotropic dissipation 
of a fluctuating velocity eddy^ cf. Tehnekes and Lumley (1974). Hence, for 
the TKE closure equation (3S) is replaced as 


V* s k'% (44) 

The dissipation length scale may be expressed in terms of the isotropic 
dissipation rate of turbulence s, cf. Hanjalic’' and Launder (1972), defined 
as 

= («) 

where is the correlation constant. Combining equations (44) -(45) yields 

V* = (46) 

which corresponds to the two-equation TKE model definition for turbulent 
effective viscosity. 

A partial differential equation system for the determination of the 
turbulent kinetic energy k and the dissipation rate of turbulence e is 
required. For turbulent incompressible boundary layer flow, the appropriate 
system is (Gebeci and Smith, 1974) 


Lfici = n— + _ -5- 


-u ® 3 k 
2k » 


- 

(.3yJ 


+ e * 0 


(47) 




V® 3e 


- + C| e2k“i = 0 (48) 
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The various are correlation constants of the model. Table i lists values 
used in this analysis, as recommended by Hanjalic'* and Launder (1972) for 
two-dimensional shear flows. 


TABLE 1 

CORRELATION COEFFICIENTS IN TKE CLOSURE MODEL 


Variable 

Equation 


V* 

(46) 

C^=0.09 

k 

(47) 

C,^=1.0 

e 

(48) 

Ce=1.3, CJ=1.44, C|=1.92 


Closure of equations (29)-(30) also requires establishment of 
appropriate boundary and initial conditions. The boundary conditions for 
solution of equations (29)-(30) are determined by inspection. At the surface 
y=0, no-slip and no injection is assumed yielding 

u(x,0) = v(x,0) = 0 (49) 

At the freestream, y>6(x), from equations (30) and (31) 

= 0 (50) 

The boundary conditions for k and s are vanishing at both the inviscid 
freestream edge and at the solid surface. For the latter location, it is 
also necessary to enforce the wall damping influence on the velocity 
fluctuations, in the manner of the Van Driest damping function oj for the MLT 
closure. 

Since equations (34), (47) and (48) also display initial -value 
character, an appropriate specification is required to initiate a solution. 
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Any arbitrary profile for u is admissable that also satisfies equations (49)- 
(50), For k and e in boundary layer type flows, cf. Launder and Spalding 
(1972), the length scale is proportional to y in the iinmediate vicinity of 
a wall. From an exact analysis, hence 

* 

Where Gq-G. 09 and ic is von Karman's constant (0.435).' Away from the wall, 
eventually becomes independent of y and levels off at a value about equal 
to Cq 6, where 6 is the local boundary layer thickness. Assuming a continuous 
distribution of between these extrema, and using the MLT model to compute 
an initial distribution of both k and e can be determined using equations 
(45), (46), and (51). 

/hi added complication of the problem specification is that 6, the 
boundary value solution domain to be spanned by the. finite element discreti- 
zation, is variable with x, see Figure 2. A transverse coordinate stretching 
transformation can efficiently compensate for boundary layer thickness growth. 
Referring to Figure 3, a useful transformation is 

y^fi(x) 

n s f ' (52) 

- fj^(x)j/f 

where n is the normalized coordinate lying between two piecewise continuous 
surfaces, f^(x) and f 2 (x)>fp that bound the solution domain R^, and where 
f is a normalizing factor. For the boundary layer flow, f^ typTcally 
corresponds to a surface and is a constant, while f 2 is appropriately 
specified. The chain rule for differentia^tion yTelds ^ „ 
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_3_ „ ^ _3_dn 

3x “ 3?"3 x 9hHx 


-L 

35 




(f2-fl)/1 


+ n 


-F"* - f 


f^-fl 


J- 

3n 


(53) 


where superscript prime denotes the ordinary derivative. The resultant form 
of the three initial -valued equations (34), (47), and (48), and the continuity 
equation (29), is 


Lm = [^-(h2.nh3)|^B 




(54) 


Uu) = 0[^ - (h2^nh3)|^]a . f 


L(k) = b[^- (h2+nh3)i^k + v-f 


^ nl 3 


3t1 



v®h 

1 ak 

3n 

Urn 

3ti 


(h2 

+ Tih3 

i®hj 

as 

- r 


'3T1 



”1ff • ■ • " 




~ 3s 
^3 t1 

3'nJ 


(56) 


(57) 


The functions h^-, l_i_3, are related to the metric of the coordinate 
transformation and defined as 


h, s 


r'l 

[rHfg-fJjhJ 


( 58 ) 


The superscript prime denotes ordinary differentiation with respect to x, and 
the h^. are at most a function of x or §. 

Finite Element Solution ATgorithiii 

The finite element solution is established assuming all •dependent 
variables and parameters are interpolated on disjoint interior subdomains 

where x is a generalized initial-value coordinate. Defining the 
numerically determined finite element approximation to the true solution as 
q*, then 

H 

q*(y»x) = I q (y,x) (59) 

e-a ® 

where M is the total number of finite element domains spanning R^. Gn each 
domain then 

qg(y.x) - qgCXj.x) 5 {\(x^)rtQ(x)}e 

where the union of S2. forms the solution domain fl s R^xx e yxx. Furtherniore, 
the elements of {N|^} are interpolation polynomials complete to degree 1* While 
the elements of {QL are the to-be-determined expansion coefficients. Using 
the familiar Gal erki n-Wei ghted Resi duals formulation, cf. Baker and Soliman 
(1977), and for the dependent variabTe q identified With u(x,y), the eTe^^^ 
form of the finite element solution algorithm for equation (34;) becomes, cf. ' 
Baker (1978)^ . , 



19 


where § Is the finite element assembly operator. In hypemiatrix form, cf. 
Baker and Solinian (1978), equation (61) for the case of no grid stretching 
becomes 

S|Ag{U}^[A3Q00HU}^ 

+ Ag[{XNUEFF}g[A-3011] + {V}gCA30Q13]{U}g 

+ AgPj{A10}] 5 {0} (62) 


In equation (62), the matrix elements of {XNUEFF} are the nodal values of the 
effective viscosity v®, equation (32), on Rg. Completion of the solution 
specification is achieved as 

{U(Xq)} u(x^j,y) (63) 

which corresponds to a mapping of the initial condition for Oi onto the nodal 

coordinates of ^u_. The rank of the global matrix system equation (62) is 
e ® 

one less than its order, to account for the no-slip boundary condition 

u(x,0) « 0. Hence, application of the finite element algorithm to equation 

(34) has yielded a system of ordinary differential equations (62) with 

initial conditions equation (63) for solution of the node point distribution 

of the discretized velocity, u, i.e. {U} = ^{U} . 

e ® 

Solution of equation (62) requires determination of the discretized 
equivalent transverse velocity distribution {V}, hence {V}_. The x dependence 
in v(x,y) is parameterized, resulting in an ordinary differential equation 
specification in the form 


3?(x,y) 

ay 



V'(x) * f(x,y) 


-{U(x)>" 


(64) 


f(x,y) 


(65) 


Hence, the y-dependence in f(Xjy) is discr’etely determined at node pdints of the 
finite element discretizatioh of R^, and equations (64)>-(65) express solution 
for transverse velocity in the standard form of an ordinary differential 
equation with X as a parameter. The trapezoidal rule rtumertcal integration 
algorithm was anployed for the present sol Utions yielding 


V(x) » V„(x) + 
n+1 " 




In equatin (66), ny - ■■ ^n^ integration step-si 2 e for equation (66) 

which lies in the domain spanned by the finite element discrettzatfon of 

equation (62). A second order accurate backwards difference formula to 
was employed to; evaluate the vector {U}'* as 

■ Sji-TolTh-^ 

p p-X' p p-i' L 

- (h„ + h„ ,)2{U}„ , 

' p p-1' P“1 

+ (h J^Wd-zI (67) 


In equation (67), hp and hp_| are the current and inunediate past Ax integratit 
step-sizes used for the sol otion of equation (62) . 

Restatement of the fi nite element al go now u|ing the coO nate , 
transformation equation (52), is obtained by simply letting pe represent each 


of Ug, kg, and Cg, in equation (60), yielding 


I crf[A3d001{Q}(' 


(^30001 + [hgfETA}^ - tAAOOOlJ; {Q} 


+ A^h,{XNUEFF}irA301il{QL + A .{SORCQL 
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The elements of the coTiunh matrix {ETAlg are the (stationary) nodal coordinates 
Of the finite element discretization of R^. In equation (68), {SORCQlg is the 
source/sink termi distinctive for each identification of q. From equation (61) 
for q = ti 


fSORGU}^ = pJtAlOl 


(69) 


which is independent of e since pj is element independent. For q = k, referring 
to equation (56) 

f v®f^]*dn «■ {XNUEFF>^ f{N^}{U}^{N|j}'{N|j}'{U}g{N|j}dn 

4 . . 4 


= ig[{XNUEFF}'’'ft406UI{U}gj{U}| 


(70) 


In equation (70), [A40bll] is a hypeVimatrix of order k+2, where k is the 
complete degree of the finite element interpolation polynomial, equation 
and each matrix element of [A40011] is similarly a square matrix of the same 
order. Hence, the source term for k can be expressed as 

{SORblOg * -|{XNOEFF}^IA4Q0111M^^^ [A200j'{EPSl^ (71) 

where the bracK have been added to emphasize that the inner matrix product 
must be performed first. 

The source tsp™^ ier* e in equation C5?^) involves the. irv?at1pnaT^^^^^f^^^ 


e/k. Since both s and k are interpolated Over R| using equation (60), and as 
an alternative to. an exact handling e/k,. an element average is employed as 


^ ^{AlOj'^IEPS? - -rr- ' . 

To wi thid a scal ar factdie> l^PRC|P>g then become s i 
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{SORCEP}g s -Cjs/kg [{XNUEFF}g[A40011]{U}gj{U}g + C|i7!T [A200]{EPS}g (73) 

The energy inner product, equation (6), for the boundary layer equation 
system is required established to evaluate convergence in the energy norm. The 
discrete solution energy norm is 

e-1 


-11 
2* Re 



{XNUEFF}g[A3Qll]j 



(74) 


Convergence properties of solely the finite element discrete solution are 
determinable as, see equations (34), (8), and (9) 

where 


nxiik ^ 




dy < 


(76) 


For the boundary layer problem, 2m is the order of equation (34), hence m=l, 
and k is the degree of the highest complete polynomial in the approximation to 
u(x,y). 

Error and convergence are also measured in terms of the boundary layer 
integral parameters. As discussed, boundary layer displacement (6*) and 
momentum (e) thickness are variables of primary interest in engineering evalua- 
tion and evaluable in terms of the p-norms. For the discrete solutions. 
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These norms, are herein yse^ pHniftriTy soTutlon accuracy in terms of 

the shape factor H, 

H = 6*/8 (79) 

and the skin friction coefficient 



Spectficany, equation (80) is evaTuated using the Ludwieg-Ti liman formuTa, 
equatidn (41). 

A Laminar Parabolic Flow 

The fundamental requirement is to assess the convergence character of 
the finite element solution algorithm applied to the non-^l inear parabolic 
eqtiatiOh (30). the elemehtary case accrues fpr 1 aminar flow* hence u^r = o. 

A further simpltflcatiOn is to assume that the transverse velocity component 
y is everWhere cOhstant, speci ff cally zero. This assumpti oh viola^^^ 
physics of the flow, except in the case u(X,y) u(y), wherein the continuity 
equation (54) is also satisfied. However, since the Tnitial-yalue splution is 
of prime interest, for the first case the continiiity equation Was discarded 
and V set to zero. For non-vanishing pressure gradient, equation (55) then 
reduces to 
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which IS the most elementary model form for the developed equation system that 
displays the essential required non-linearity. The boundary conditions for 
solution of equation (81) remain expressed by equations (49)-(50), where 5 is 
now assumed a symmetry plane, and the initial condition is the slug profile 
illustrated in Figure 4. Following an extensive numerical test program 
(Soliman, 1978), this initialization for u was found mandatory to eliminate 
the initial data as the primary error source in the energy error norm, equation 
(74). Selecting the distribution illustrated provided a uniform initial energy 
for all k, l5k< 3, on both uniform and non-uniform discretizations of A, the 
span of R^. The particular data set selected was U^= 300 ft/sec, Re = 0.7x10®/ 
ft and Ap = 1.458 X 10"®. Following experimentation, quantization of convergence 
was determined facilitated after marching the solution downstream for an axial- 
distance of approximately 0.5 ft; therefore, all data presented were measured 
at Ax = 0.8 ft. 

Figure 5 presents computed solution error in the energy norm as a 
function of uniform discretization refinement for linear, quadratic, and cubic 
finite element interpolations. The numerical results confirm the prediction of 
the linear theory, equation (74), for jc = l and 2. Specifically, convergence 
is exactly quadratic for the linear finite element solution, and essentially 
fourth-order for the quadratic elements. In contrast to theory, convergence 
is also fourth-order for the solutions obtained with cubic finite elements. 

Note that the level of error capable of quantization is of the order of 1Q“®, 
which has been confirmed as the upper bound on accuracy for the equation 
solved, which is non-linear. As will be documented throughout this test 
program, the coarse grid linear finite element solutions display accuracy 
superior to that predicted by strict adherence to the convergence curve. 
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Conversely, those solutions obtained with quadratic finite elements display 
coarse grid inaccuracy. The finest grid case reported for a quadratic finite 
element soTuti on (H = 40) also displays a somewhat larger error than predicted 
by the fourth-order convergence rate, which is further confirmation of the 
accuracy level attainable in actual practice. As noted in the header of 
Figure 5, convergence is monotonic and from above for the linear and cubic 
element solutions, while the quadratic element solution displays oscillatory 
convergence. 

Figures 6-7 present computed accuracy and convergence in the engineer- 
ing norms of shape factor and skin friction, equations (7/)-(80), for k = l, 2 
and 3. Convergence is quadratic for the linear finite element solution, and 
from above in the case of the shape factor and from below for skin friction. 
Fourth-order accuracy is displayed by the quadratic finite elements, and in 
contrast with the energy norm, the coarse grid solution demonstrates a super- 
accuracy. The slope of a straight line drawn between the two data points for 
the cubic element solution is 5.8, which is in close agreement to the six that 
is predicted by the linear theory, equation (7), or specifically k = 3 and m = l 
in equation (75). The next consistent cubic element discretization requires 
M=54, the results for which far exceed the quantizable error of 10“® for this 
problem. Hence, the cubic element formulation becomes essentially impractical 
for equations of this type in fluid mechanics. 

The influence of employing a non-uniform finite element discretization 
of Rl", on the computed error in the energy norm, was determined and results 
are presjOnted in Figure 8 for k = I and 2. The abscissa for this curve is now 
the largest finite element on Rl. The non-uniform discretization increases 
by up to a factor of two: the level of absolute error in the energy norm for the 
iTnear element solution, k = l. Alternatively, for the quadratic, the error 
"level is less than that associated with the uniform discretization, with an 
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indication that an optimum grid exists for which an essential minimum exists 
for energy. Figures 9-10 provide the comparison on the basis of error in H and 
Cf. In these nonns, a favorable effect is accorded use of non-uniform 
discretization. The absolute error is substantially decreased using a non- 
uni fom discretization for k = l, and there is an optimum grid associated with 
an absolute error minimum in both shape factor and skin friction. The levels 
o.f absolute error are uniformly decreased for the k = 2 solutions, as well, 
although no optimum discretization is evident. The non-uniform discretization 
results display an essential fourth-order accuracy in shape factor norm, while 
a nominal third-order accuracy is evidenced for skin friction. In all cases 
then, the use of a non-uniform discretization is not contraindicated, and subse' 
quent results for a physically neaningful equation statement will confirm this 
indication. 


Laminar Boundary layer Flow 

The second and physically meaningful test case corresponds to incompres- 
sible laminar boundary layer flow impinging on a sharp leading edge with 
pressure gradient. The nominal freestream velocity remains I 4 , = 300 ft/sec and 
Re = 0.7xl0®/ft. The slug-profile for u. Figure 4, remains initial condition 
for the stream-wise velocity, while the transverse velocity v is assumed zero 
at solution initiation. The transverse velocity remains zero until sufficient 
solution information has been generated to facilitate evaluation of the backwards 
difference formula, equation (66). To maintain a uniform evaluation of the 
initial data energy norm, equation (74), for all k and all initial-value matrix 
assembly operators the node at the knee of the initial slug profile for the 
streamwise velocity is kept at the same position (0.2A) for all grid refinements. 
To facilitate the required comparison with a popular finite difference solution 
algorithm, numerical results are also obtained using the equivalent of the 
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familiar Crank-Nicholson algorithm for equation (62). This is achieved by 
rendering diagonal the initial value matrix [A3000], within the linear finite 
element formulation. The specific matrix equivalent for the initial -value 
matrix [A3000]g associated with Crank-Nicholson is 



For illustration, upon completing the assembly operation, defined as S(^ s 3 
herein, equation (62) can be reexpressed on a uniform grid in the (finite 
difference) recursion form 


“j“5 * fe+i ■ “j-l] 




(83) 


Upon application of the trapezoidal integration formula for u^, the resultant 
algebraic equation system is identical with the Crank-Nicholson algorithm, 
cf. Roache (1972). 

Accuracy and convergence evaluations were again obtained at an axial 
displacement of 0.8 ft downstream from the leading edge of the plate. For 
linear finite element functions, k = l equation (60), convergence with discre- 
tization refinement in the energy norm is computed uniformly quadratic, for 
Initial -value matrix structures S = 0, 2, and 3, see Figure 11, and with 
negligible data scatter. Herein, S = 0 corresponds to the exact finite element 


formulation, while S-2 corresponds to the diagonalized structure previously 
employed by Baker and Manhardt (1977, 1978). (The operation S=1 is another 
diagonalizing operator that is consistent only for linear element formuTationS:^ 
and thus reli gated to history.) The errors are calculated with respect^ 
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estimated exact value of the energy norm, based on the assumption that the S = 0 
fine grid solutions do indeed converge quadratically. The lowest error in the 
energy norm is accorded the finite element solution algorithm, S = 0, and the 
convergence is from below. Convergence is from above for the alternative 
diagonalized formulations S = 2,3. This confirms prediction of extension of 
the theoretical convergence theorem, equation (16), to this practical non-linear 
problem, and the finite element solution indeed minimizes the error in the 
energy. Note also that the error for the coarse grid finite element solution 
(M=10) falls below strict adherence to the convergence curve, confirming and 
firmly quantizing existence of coarse grid accuracy as previously predicted by 
Popinski and Baker (1976). Convergence in the engineering norms of K and Cf 
is also firmly quadratic for all S = 0, 2, and 3, see Figures 12-13, with the 
finite element solution over-predicting shape factor and under-predicting skin 
friction. Tliis trend is again reversed for S = 2 and 3 solutions. The 
difference in level of error is much less pronounced in these norms, although 
the S=0 error is uniformly minimum, and coarse grid accuracy remains evident 
in the engineering norms. 

Accuracy and convergence evaluation for linear, quadratic and cubic 
elements, l^k<3, for the same test case, and as measured in the energy norm 
for the consistent finite element initial -value matrix S^O, are presented in 
Figure Id^. The solid curves are of nearest integer slope, and the demonstrated 
conver*gence rate for k=l, and the coarser grid solutions for k = 2, predicts 
extension of the linear theory for both k to this non-linear equation. The 
cubic finite element results, k = 3, and the finer grid solutions for the 
quadratic functions, k = 2, fail to adhere to the convergence curves for error 
in the solution less than about lO’’^. This further confirms existence of a 
practical bound on actual perfonnance of high order accurate numerical solution 
algorithms for non-linear equations of the boundary layer type. Furthemore, 
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the performance of the cubic element formulation is truly dismal. By compari- 
son, absolute numerical accuracy for solution of linear parabolic equations was 
reported better than 10"^^, see Baker and Soliman (1978). Only for linear 
elements is convergence monotonic and from below. For the k = 2 and k = 3 cases, 
convergence is oscillatory and starts from below. 

As a measure of solution economy. Figure 15 presents the error in the 
energy norm as a function of discretization refinement now expressed as the 
number of nodes on the solution domain. The abscissa is equivalently the rank 
of the Jacobian associated with solution of equation (62), and represents the 
amount of computational work required to obtain the solution. For all cases 
tested, except perhaps for the coarsest grid, the quadratic element solutions 
demonstrate a definitely superior economy and accuracy in the energy norm. The 
cubic element formulation is even less favorable on this comparison basis. 

Figures 16-17 present computed convergence in the engineering norms for 
the finite element solutions obtained for k = 1, 2, and 3. Except for the 
coarse grid solutions obtained using the quadratic elements, convergence is 
generally quadratic for all k. The exception is the coarse grid solutions 
obtained with quadratic elements, wherein the M=10 solution is super-accurate, 
exceeding even 4^^ order convergence as indicated. This is most probably the 
direct consequence of the convergence being oscillatory. Note that the 
degraded convergence performance for 2£k<3 occurs at error levels less than 
10“5. As was determined in the energy norm comparisons, a practical upper 
bound on attainable solution accuracy definitely exists for this non-linear 
parabolic equation, and appears to be of the order lO"^. The overall superi- 
ority of the quadratic element formulation remains apparent, with the associated 
ei ror roughly an order of magnitude less than that of the linear element 
solutions for approximately the same amount of computational work (same number 
of nodes). 
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These laminar flow results were obtained using unifom finite element 
discretizations of R^. It is necessary to investigate the important effects of 
a non-uniform discretization, which is a key feature required for efficient and 
accurate turbulent flow computations. A smooth progression of non-uniformity 
of the finite element discretization is desired, and is attainable using a 
geometric progression to locate nodal coordinates on R^ as 


Vl ° TI 


Ri (e-l) 


l<e<M 


1 


(84) 


In equation (84), y^^j^ is the extremum nodal coordinate of R^, y^ is the 
coordinate of the first node of R^, p is the geometric progression ratio, and M 
is the total number of finite elements R^ spanning the domain R^. 

The M=80 linear element discretization was chosen as the base case 
for comparison of all non-uniform discretization results. The node at the knee 
of the strearawise velocity slug initial profile was maintained at a constant 
coordinate (20% of the domain), such that the initial data energy norm 
evaluation was uniformly a constant for all non-uniform discretizations. Figure 
18 present computed accuracy and convergence in the energy norm for the linear 
finite element algorithm for a range of pressure gradients that yielded a zero 
and a ±50% change in solution energy compared to the initial data evaluation. 

The non-uniform discretization solutions are represented by partially shaded 
symbols, and the sign next to the symbol indicates that the sign of the norm 
changed with respect to the estimated correct value in comparison to the uniform 
discretization results. The computational advantage of using a non-uniform 
discretization is clearly demonstrated; the error in the energy norm is minimum 
for any number of elements spanning the solution domain using any non-uniform 
discretization. Furthermore, the non-uniform convergence curve, shown as a 
dashed line, passes through a minimum (zero) which indicates existence of an 
optimum grid (Ha»27, p**1.2) for this particular problem. 
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In the engineering norms, the effect of pressure gradient on the error 
Is more apparent. Referring to Figures 19-20, the Increase In level of error In 
the shape factor and skin friction is approximately five for the intermediate 
pressure gradient and ten for the extremal pressure gradient plotted. The 
beneficial effect of use of a non-uniform discretrization on error in the 
engineering norms decreases as the pressure gradient increases, and the discre- 
tization which extremizes the energy does not necessarily yield the lowest error 
in the engineering norms. However, in all cases, the accuracy attainable using 
a non-uniform grid is universally superior in any norm. A non-uniform grid 
containing approximately 25 elements will always yield accuracy comparable or 
superior to the 80 element uniform grid case at a factor of 3 reduction in 
computer CPU. Since the data collapse to an essential single curve. Figure 18, 
the energy norm appears the superior mathematical measure of accuracy for this 
non-linear parabolic equation system. 

The comparison is required established for the alternative initial-value 
matrix operators, S = 2 and 3. Figure 21 presents computed convergence in the 
energy norm for the selected range of pressure gradients for solutions obtained 
employing the diagonalized initial -value matrix S = 2. Contrary to the numerical 
experience with the consistent form S = 0, see Figure 18, the absolute error in 
the energy norm decreases with the Increase In pressure gradient level. The 
convergence rate remains essentially quadratic for uniform discretization 
solutions, but there is select data scatter and some evidenced coarse grid 
inaccuracy. Only in the case of zero pressure gradient is the error in the 
energy norm consistently reduced using a non-uniform discretization. The abso- 
lute level of improvement is drastically reduced in comparison with the consist- 
ent assembly results. The use of a non-uniform discretization modestly increases 
the error in the energy norm for the intermediate pressure gradient, Ap = 1.525 
It 10“S, while converging from above. A select modest improvement, coupled with 
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convergence from below, is accorded the extreme pressure gradient solution, 

Ap = 3.056x 10“®, for non-uniform grid comparison. 

Convergence with discretization refinement, in the engineering norms 
for the diagonal matrix S = 2, is presented in Figures 19-20. Quadratic converg- 
ence is confirmed in all cases for the uniform discretization results, with 
negligible data scatter, and the absolute error in the norm now increases with 
the increase in pressure gradient, in accord with the results obtained for S = 0. 
In distinction, however, note that use of a non-uniform grid for zero-pressure 
gradient exerts no consequential effect on solution accuracy in either engineer- 
ing norm or in energy. In contrast with error in the energy norm. Figures 21, 
non-uniform discretization can reduce error in both the shape factor and skin 
friction for non-zero pressure gradients. The absolute level of error in all 
three norms, comparing the S-0 and S*2 results, is approximately comparable 
to the uniform grid results. However, the S«0 results are clearly superior in 
all norms for the -ase of zero pressure gradient. 

The accuracy and convergence performance of the Crank-Nicholson finite 
difference initial-value equivalent, S = 3, is not consequently distinct from 
the S=2 results. Figures 24-26 present the corresponding computed error in 
the solution norms as a function of the discretization refinement. Convergence 
is again essentially quadratic, with only modest data scatter for uniform 
discretization, with the absolute error higher than that associated with the 
finite element solution (S = 0) for zero pressure gradient and almost the same 
for non-zero pressure gradients. Somewhat improved accuracy accrues to use of 
non-uniform discretizations, in contrast with the diagonalized matrix (S = 2) 
results for zero pressure gradient. Non-zero pressure gradient performance is 
nominally identical. In terms of error in shape factor and skin friction, the 
S“3 and S-2 results do not differ. Interestingly, the use of non-^uniform 
discretizations for zero pressure gradient again does not improve solution 
accuracy in either engineering norm. 
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Based upon these results, the consistent finite element initial-value 
matrix (S = Q) form for the linear (k=l) solution algorithm demonstrates 
consistently superior solution behavior in terms of accuracy, convergence and 
economy for this practical non-linear parabolic equation system. Previous 
results (Baker and Soliman, 1978) indicate this to hold as well for use of higher 
degree (k>l) finite element polynomials for a linear parabolic equation. Hence, 
numerical evaluation of accuracy and convergence obtained using higher-degree 
finite element functions is conducted for S = 0 only. Shown in Figure 27 is the 
computed solution error in the energy norm as a function of discretization 
refinement. The fourth order accuracy of the algorithm, which was documented 
earlier for zero pressure gradient (Figure 14), is now degraded to second-order 
for the cases involving non-zero pressure gradient. A plausible explanation for 
this is that the oscillarory convergence in the energy norm, associated with the 
zero pressure gradient solutions, changed to monotonic convergence from below 
upon applying a pressure gradient. Furthermore, coarse grid inaccuracy is now 
evidenced for non-zero pressure gradients, since the corresponding data points 
lie above the convergence curve. The absolute level of error for the non-zero 
pressure gradient case is five times smaller then that associated with the 
linear finite element solution for the finer grids (see Figure 18). This 
improvement over the linear element solution degenerates, however, as the grid 
progresses to coarse. The error for the M=5 quadratic element solution is 
basically Identical to that for the M=10 linear element grid. Use of a non- 
uniform discretization, selecting the M=40 element grid as the base case, 
consequentially reduces the absolute error level for all pressure gradient 
including zero, in agreement vyith the k-1 solutions. In clear distinction, 
however, an optimum grid that extremizes the energy cannot be detected. 

Figures 28-29 present the corresponding data on quadratic element 
solution error measured in the engineering norms. Convergence is essentially 


(^uadrittc; faf unifofnl dtscfetfzattort^^ Witfi cbaPse grids display fourth-order 
convergence for all pressure gradients* OsciTl atory convergence Is the general 
trend for all the cases, except for skin friction v/lth zero pressure gradient. 
The quadratic finite element solution is relatively favorable for the cases with 
pressure gradient, since the error in the shape factor is two orders of magni- 
tude smaller than the solution error obtained with linear finite elements. This 
compares to only one order of magnitude difference for the zerti pressure 
gradient case. Non-uniform discretizations display improved absolute error 
level for all cases, and there is no. indication of an optimum grid. These 
solution convergence trends are unchanged when measured in the skin friction 
norm, except that the error now decreases with an increase of pressure gradient. 
Use of non-uniform discretization again decreases the absolute level of error 
with no indication of an optimum grid. 


A Turbulent Boundary Layer Flow 

Acceptable resolution of near wall damping phenomena is an essential 
key feature of turbulent flow computations. Since use of a uniform discreti- 
zation would require an impractically large number of elements to span the 
solution domain, a non-uniform finite element discretization is required in 
all instances to obtain satisfactory computational efficiency in concert with 
acceptable solution accuracy. Baker and Hanhardt (1977b) have determined 
that linear element solution speed and accuracy, using an explicit integration 
algorithm and S = 2, both accrue using a finite element discretization with 
nodal coordinates determined according to the geometric progression equation 
(84). Based upon the laminar flow results discussed, and assuming the linear 
equation theory extensible to the more non-linear turbulent flow equations, 


the linear finite element algorithm should yield a quadratically-convergent 


procedure. Similarly, by extension, the use of quadratic elements should 


yield 4 fourth-order convergent algorithm in energy. 
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The first requirement in this analysis is to confirm Indeed that the 
developed finite element algorithm is capable of accurate prediction of 
turbulent flow for which comparison results exist. This is provided solely 
by experimental data, and a particularly challenging configuration corresponds 
to the IDENT 2400 data, reported in the proceedings of the AFOSR-IFP -Stanford 
Conference on Computations of Turbulent Boundary Layers (1968). IDENT 2400 
is the Bradshaw relaxing flow data set, which corresponds to evolution of a 
non-equilibrium subsonic boundary layer induced by abrupt removal of a 
moderately adverse pressure gradient from an initially equilibrium flow. 

Nominal freestream velocity (Um) is 33.5 m/s, wind tunnel background 
turbulence level was less than 0.155, and the reference unit Reynolds number 
is 2.38x 107 m“i. The test case is considerably demanding since non-equilibrium 
phenomena are involved in the relaxation process. The base case results were 
generated using the linear element (k = l) algorithm and a non-uniform discreti- 
zation. Following considerable numerical experimentation, an adequate 
resolution of the wall region damping was determined captured using H=30 
linear elements spanning approximately 1.56, and a geometric progression ratio 
of p = 1.222, see equation (84). Turbulence closure for the base case was 
accomplished using mixing length theory (MLT), with the parameters tc and X 
equated to their standard values of 0.435 and 0.09 respectively. For boundary 
conditions, both u and v vanish identically at the plate surface, and 3u/8y 
vanishes for y > 6, The first member of IDENT 2400 data set was interpolated 
at the nodes of ^R^ to generate the initial distribution for u, and v was 
assumed zero until sufficient data was generated to initialize the continuity 
equation solution, see equations (64) -(67). Shown in Figure 30 are comparisons 
between data and the computed solutions, for the important boundary layer 
parameters, and as obtained using the three initial -value matrix structures, 

S = 0, 2, and 3. The computed results were matched with the data at the 


second experimental profile,: as reconpeiided in the Proceed These solutions 
Were generated from equations (54)-(55)^^^^u the transfomed coordinate system 
with 20^ grid growth over the solution range of 1,3 m. Agreement with data is 
generally good, indicating the basic algorithm capable of accurate resolution 
of the physical problem. The correct trends and local extrema in 6* are 
predicted; however, the overall level of the solution curve is somewhat high. 

The level of the curve for the standard finite element structure, S = 0, is 
closer to the data than that predicted by S = 2 and 3, The computed extremum of 
v®/v, equation (32), for this case was 900, which indicates a high level of 
turbulence. 

Shown in Figure 31 is a comparison of the computed energy norms for S = 0, 
2, and 3. Note that the energy norm is minimized by the finite element solution 
S = 0 tnroughout the solution range, which generally predicts extension of the 
linear theory, equation (16), to this highly non-linear problem class. Figure 
32 presents comparison between select computed velocity profiles and data at 
three downstream stations, and agreement is generally excellent. 

To investigate the influence of discretization refinement on solution 
accuracy, the nimiber of linear elements was doubled to M = 60 while retaining 
the first node off the wall at the same physical location, to preserve 
satisfactory resolution of near-wall damping. The resulting progression ratio 
p for this non-uniform discretization was 1.089. Shown in Figure 33 are 
comparisons between data and computed solutions using the standard H = 30 linear 
element discretization and the M = 60 element discretization for the two 
different ratios of grid nodal progression. There is essential overall agree- 
ment between the two solutions obtained using 20% grid growth. The M = 60 
element discretization with 50% grid growth is in slightly better agreement with 
data at the further downstream stations except for a tendency to over-predict 
0 and 6*. The computed energy norms for these three solutions are presented in 
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Figure 34. Note that the energy is minimized by the M = 60 element discretiza- 
tion in the near field part of the curve. Thereafter, the knee in the curve 
is associated primarily with the inclusion of an extra element in the solution 
domain due to the boundary layer growth. Using S0% grid growth to keep the 
boundary layer edge within the same element throughout the solution range 
resulted in the higher overall level for the energy norm associated with 
correspondingly larger element spans. 

Figure 35 presents comparison between computed solutions obtained using 
the H = 60 linear element discretization and a M = 30 quadratic element 
discretization. The first node off the wall was maintained at the same 
physical location for both cases, which results in a progression ratio of 
1.188 for the M = 30 quadratic element discretization. Note that using a non- 
uniform discretization for the quadratic element case results in placing the 
element vertex nodes in the geometric progression while the interior nodes 
remain located at the raid-span of an element. The grid growth was 20% for the 
linear and 30% for the quadratic element solutions. In comparison with the 
experimental data, the solution using quadratic finite elements yields 
generally more accurate estimates for the boundary layer parameters than that 
obtained with linear finite elements. This cannot be directly confirmed, 
however, from noting the results presented in Figure 36. The energy norm 
calculated using quadratic finite elements has a higher level than that calcu- 
lated using linear finite elements throughout the solution range except for a 
small portion at the beginning. This is in part a direct result of a higher 
estimate of the boundary layer thickness 6, for the quadratic element solution, 
which yields correspondingly higher values of the effective viscosity v®/v. 

The computed extremum of v®/v for quadratic element solution was 937 compared 
to 891 for the linear finite element solution. 


All result? discussi^^ a fixed integraMprj step size 

X = bi05 ft, with reevaluetipn of the Jacob tWelvd integp steps. 
Thia solatibhs; w6t% as 5 umed::cQ[n\ferg^^^ dependent 

ble (streamwise Af6l^^ / u/Upa^ was unifbnnlY less thap the cbnyergence criten 

Using: a as small as lQf”f did hot attef the 
signification digit in the solution normsi Table 2 summarizes the results of 
numerical experimehts carried out to aSsess the efficiency of the algorithm for 
the Bradshaw relaxing flow test case as obtained with the M“ 30 linear element 
non-uniforin discretization with progressioh ratio p- 1.222. The tabulated 
results correspond to the final solutions at Ax = 1.3 m. The results tabulated 
for the energy nom, shape factor and skin friction show the significant place 
of the integration truncation error, as confirmed by a higher-order accurate 
solution obtained using Richardson extrapolation. The slash isolates the 
significant digit in each norm, with the upper result corresponding to the 
more accurate one obtained using half the regular integration step size (the 
Richardson step). 

The reference solution in this comparison (case 1) was obtained using a 


fixed uniform integration step size Ax = 0.05 ft with reevaluation of the 
Jacobian every twelve integration steps, which required evaluation of the 
Jacobian 66 times throughout the solution range. Using twice the integration 
step size and reevaluating the Jacobian every 12 steps, case 2, yielded 
identical values for the norms while reducing the number of passes and 
accordingly the CPU by 335J. The effect of utilizing the coordinate transforma- 
tion equation is documented by case 3, wherein the solution domain was allowed 
to grow linearly in the streamwise direction in such a fashion that the span 
of the solution dGmain at the final Integration station was 2035 larger than at 


solution initiation. This solution minimized the energy norm while the 


computed difference in the 


factor is 0.6% and in skin friction is 2.555. 


TABLE 2 

EFFICIENCY OF THE LINEAR FINITE ELEMENT ALGORITHM - BRADSHAW RELAXING FLOW (MLT) 


Grid 
Case Growth 
No. % 

Integration 
Step Size 
Ax 
ft 

Convergence 

Criteria 

e 

CPU* 

Number of 
Jacobian 
Number of Reevalua- 
Passes tions 

Percent 

Increase 

in 

Step Size 

E 

(10-3) 

H 

Cf 

(10-*) 

1 

0 

.005 

10“6 

1,00 

2341 

66 

0 

.572/35 

44 

1.389/39 

40 

.1784/91 

80 

2 

0 

.010 

10" 6 

.67 

1704 

33 

0 

.572/05 

13 

1.3896/49 

58 

.17836/48 

26 

3 

20 

.005 

10"6 

1.42 

2418 

66 

0 

.529/52 

60 

1.3805/04 

15 

.18293/30 

51 

4 

0 

.005 

10-6 

1.03 

2473 

0 

0 

.572/36 

47 

1.3894/02 

11 

.17848/87 

56 

5 

0 

.005 

10‘^ 

.53 

905 

66 

0 

.572/31 

40 

1.3893/54 

66 

.17850/20 

00 

6 

20 

.005 

10" 6 

.92 

1438 

27 

10 

.517/37 

40 

1.3721/86 

90 

. 18421/68 
70 


*Normal1zed on case number 1. 


Case 4 is identical to case 1 except that the ini ttal 
retained throughout the solution range. The difference in the energy norni 
and both engineering norms between these two cases, is beyond the significant 
digit based on Richardson extrapolation. Retaining the initial Jacobian 
resulted, however, in a 3% increase in the number of passes and GPU. 

The influence of a relaxed convergence criteria e is documented in 
case 5. Reducing e by two orders of magnitude to resulted in reducing 
the number of passes by 6Q% and a 47% saving in CPU. Kith this favorable 
economy feature, the change in the energy norm and the engineering norms 
from the reference case is again beyond the acceptable significant digit. 

In test case 6, the integration time step Ax was increased by 10% every 
time the Jacobian was reevaluated. This procedure reduced the number of 
passes by 38% and the number of Jacobian reevaluations by 60%. The energy 
norm was minimized while the change in the shape factor was 1% and that in 
the skin friction was 3%. 

The results of an assessment of accuracy and convergence trends for 
linear element solutions are presented in Table 3. These results were 
obtained employing the finite element matrix S = 0 with a convergence criteria 
e of 10“®. The span of the first element was 0.27 xio-^ ft, for the M = 30 
element and the first M = 60 element discretization, while for the second 
M = 60 element and the M=120 element discretizations, Aj^ was 0.21 xlO"^ ft. 

The larger negative value for the change in the energy. norm, normalized by 
the initial energy (AE/E), indicates a greater minimization of the energy, 
since the energy norm decreases as the solution is marched downstream. On 
this basis, solution accuracy Increases with discretization refinemeht. Note 
also that the normalized change in the shape factor (AH/H) is not affected y 
by discretization refinanent, and thus could not he used to assess conver^^^ 
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TABLE 3 

ACCURACY AND CONVERGENCE OF THE LINEAR FINITE ELEMENT ALGORITHM - BRADSHAW RELAXING FLOW (HLT) 


Number 

of 

Elements 

Progression 

Ratio 



CPU* 

Number 

of 

Passes 

AE/E 

AH/H 

30 

1.222 

.0054 

.20 

1.00 

635 

-.167 

-.036 

60 

1.089 

.013 

.085 

1.98 

638 

-.201 

-.036 

60 

1.095 

.094 

.094 

2.02 

S55 

-.204 

-.038 

120 

1,039 

.039 

.039 

4.74 

847 

-.219 

-.038 


*Normal1zed on 30 linear element grida 
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Sunnnarlzed In Table 4 Is the corresponding assessment of accuracy and 
convergence trends for the quadratic element solutions for the same test 
case. The span of the first element was chosen to be twice that of the 
corresponding linear element discretization, with twice the toal number of 
elements in the solution domain. This procedure resulted in placing the 
first node off the wall in the same physical location for the linear and 
quadratic element cases, in an attempt to maintain consistent resolution of 
the near-wall damping. The CPU time is approximately the same for the 
corresponding linear and quadratic discretizations. The solutions employing 
quadratic elements display convergence in the energy norm with discretization 
refinement, as evidenced fay a superior minimization of the energy norm. As 
for the linear element solutions, the normalized change in the shape factor 
was not affected by discretization. Comparing results in Tables 3-4 shows 
that the M=15 and M = 60 quadratic element discretizations yield a superior 
energy niinlmizatlon than the corresponding M=30 and 120 linear element 
discretizations. This is not valid, however, when comparing the M=30 
quadratic element discretization to the M = 60 linear element discretization 
results. The influence of the progression ratio used to define the non- 
uniform discretization, on finite element solution accuracy, is shown in 
Table 5. The progression ratios which yield the largest negative value of 
AE/E, i.e., extremum minimization of the energy norm, are also those which 
required the least number of passes, fhe computed effective viscosity at 
the first node off the wall (v®/v) increases as the span of the first 
element Aj^ increases, and the best results were obtained when v®/v was 
approximately equal to 2. The normalized change in the shape factor 
decreases monotonically as the progression ratio decreases; hence, it could 
not be used to indicate the preference of any progression ratio over the 
others. 


TABLE 4 


ACCURACY AND CONVERGENCE OF THE QUADRATIC FINITE ELEMENT ALGORITHM - BRADSHAW RELAXING FLOW (MLT) 


Number 

of 

Elements 

Progression 

Ratio 

Aj/Ag. 

^maX/5 

CPU* 

Number 

of 

Passes 

.AE/E 

AH/H 

15 

1.506 

.0049 

.45 

1.00 

608 

-.178 

-.038 

30 

« 

1.188 

.014 

.16 

1.96 

614 

1 

• 

H-a 

CD 

-.038 

30 

1.200 

.0087 

.20 

2.01 

638 

-.193 

-.037 

60 

1.079 

.023 

.076 

4.85 

649 

-.229 

^•038 


^Normalized on 30 linear element grid (Table 3). 
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TABLE 5 


INFLUENCE OF PROGRESSION RATIO ON LINEAR FINITE ELEMENT SOLUTION ACCURACY - 

BRADSHAW RELAXING FLOW (MLT) 



Progression 

Ratio 



NmW)er 

of 

Passes 

v®/v 

AE/E 

AH/H 


1.222 

.0054 

.20 

635 

1.151 

-.167 

-.036 

1.211 

.0069 

.20 

. 617 

1.313 

-.183 

-■034 

1.200 

.0088 

.20 

' 605 

.1629 

-.220 

o 

* 

1 

1.189 • 

.013 

.16 

606 

2.205 

o 

Csl 

Csi 

1 

- .033 

1.178 

.017 

.16 

610 

3,181 

-.163 

1 

« 

o 

ro 

1.167 

.021 

.15 

. 623 

4.799 

-.163 

-.020 

1.155 

,027 

.16 

646 

7.387 

-.187 

-.014 
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Accuracy and Convergence for Turbulent Boundary Layer Flow if 

i-i 

The presented results document viability of the finite element algorithm Dj 
and the discretization philosophy for acceptably accurate turbulent boundary 

V 

layer flow prediction. A tightly controlled numerical test case, analogous to 
that employed for the laminar flow analysis, is required to quantize accuracy 
level and convergence with discretization refinement. The case corresponds 
essentially to transition to turbulent flow of the laminar Slug start in zero 
pressure gradient. The test conditions were selected identical to the Wieghardt 
data set (IDENT 1400, Proceedings of the Stanford Conference (1968)) with 
constant freestream velocity (l^ = Uj = 33 m/s) yielding a unit Reynolds number 
of 2.19 X 10® per meter. Five different non-uniform discretizations were used to 

V-' 

study accuracy and convergence with discretization refinement. The total number I; 
of elements H spanning the solution domain and the corresponding node 
progression ratios p are listed in Table 6 for the linear and quadratic finite 
element solutions. All computed solutions were initialized essentially 
identical to the experiment, wherein a turbulence-free uniform flow impinged 
upon the plate leading edge, using the slug start profile shown in Figure 4. 


TABLE 6 

DISCRETIZATION DATA - TURBULENT FLAT PLATE FLOW 








H 

Ir - 1 

12 

24 

36 

48 

60 

“ 1 

P 

1.627 

1.222 

1.125 

1.083 

1.061 

M 

6 

12 

18 

24 

30 


k«2 


p 2.8,14 1.510 1.271 1.176 1.110 
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As before, the first node off the plate was held at the same physical position 
for all discretizations. The number of elements between the plate and the 
knee of the velocity profile was always one-sixth the total number of elements 
spanning the solution domain. No turbulence transition model was employed; 
instead computational transition from laminar to turbulent flow was specified 
to occur when shape factor H achieved 90^ of the fully developed laminar flow 
value. 

Figure 37 summarizes computed solution error obtained with the linear 
element algorithm as a function of discretization refinement. Convergence In 
the energy norm Is essentially quadratic for the three Initial -value matrix 
structures, S = 0, 2, and 3, with the finite element algorithm S = 0 again 
yielding the smallest level of error for any M. Note in all cases that 
convergence Is from above. As in the laminar flow results, the finite element 
algorithm displays accuracy for the coarse grid that Is superior to strict 
adherence to the convergence curve. Convergence in shape factor. Figure 38, 

Is essentially quadratic and from below for the diagonalized (S = 2) and Crank- 
Nicholson algorithm (S = 3). The finite element results do not display a 
convergence trend In shape factor. Specifically, the error is uniformly 
constant and smaller than that for either the S = 2 or 3 results. Since shape 
factor Is the ratio of 6* and 0, see equations (77) -(79), their convergence 
properties were measured. As shown In Figure 39, convergence in both S* and 
0 Is quadratic on coarse grids and nearly fourth order for finer discretiza- 
tions. Since the curves are parallel, the error In the shape factor remains 
essentially constant as determined. Figure 40 shows that convergence measured 
In skin friction is essentially quadratic for all three forms S - 0, 2, and 3, 
with the finite element solutioi yielding the smallest error level for any M. 

Based upon the experience with laminar flows, a fourth order accurate 
algorithm Is anticipated to result from use of quadratic finite elements 
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assumtng the Tinear thapry holds. It 1s of parttcuTajf^ irtiportariee^^^t^^ ascertain 
this, since the energy norm is now complexly related to level of turbulence 
within the flow (evolution) through the eddy viscosity. Recalling the energy 
definition, equation (6), and that eddy viscosity involves u shear, the 
specific form of the energy norm is 



From the Bradshaw test results, recall that v®/v ranged to 10^. Since the 
effective diffusioh cbefficieht itself is strictly dependent Ujjon the computed 
evolution of u, the non-linearity of the subject equation system will exert a 
profound impact on the convergence evaluations. 

Figure 41 presents the error computed in the energy norm as a function 
of discretization refinement for the linear and quadratic element algorithms. 
Convergence is from below and of generally fourth degree for the quadratic 
element solution, which predicts extension of the linear theory. However, 
the results for the finest discretization show a significantly larger absolute 
error than predicted by the convergence curve. This is interpreted again as 
an indication of the limit of practically useful discretizations. The accuracy 
of the quadratic element solutions can be a factor of up to 50 improvement 
over the corresponding linear element results. Figures 42-43 present error in 
the engineering norms as a function of discretization refinement. Convergence 
is oscillatory in both norms for quadratic elements, as experienced in the 
case of laminar flow, and of essentially fourth degree to the attainable limit 
of accuracy* The absolute error in the engineering norms for the finer 
discretizations is considerably larger than predicted by the convergence curve, 
confirming the experience in the energy norm. 
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Accuracy Evaluation Using the TKE Closure 

As a summary computational study, the turbulent kinetic energy two- 
equation closure model was evaluated using the finite element algorithm, with 
primary emphasis on solution economy. The test case corresponds to the 
Bradshaw data set discussed previously. Details on solution initiation are 
given by Soliman (1978), and consistent accuracy and convergence trends were 
computed using the TKE closure model for both the linear and quadratic element 
algorithms. 

The efficiency of the solution algorithm employing the TKE closure 
model can be appreciably improved by using one Jacobian for the three depend- 
ent variables, resulting in a considerable reduc'd on in required memory 
storage. Table 7 summarizes comparisons between different methods of handling 
the Jacobian- The reference solution (case 1) was obtained using the correct 
Jacobian for each of the three dependent variables u, k, and e. Employing the 
u Jacobian for each dependent variable solution resulted in deterioration of 
accuracy, as evidenced by the larger value of AE/E, and an overall B% increase 
in CPU. The third solution was obtained using the k Jacobian for each of the 
three dependent variables. This shows an improvement in accuracy over the 
three Jacobian reference case, as evidenced by a minimum AE/E and a 9% saving 
in computer CPU. No specific trends were indicated in AH/H. To investigate 
influence on solution accuracy, of the accuracy of the turbulent viscosity • 
evaluations within the Jacobian, was deliberately under-evaluated and 
convergence of the matrix iteration evaluated. Case 4 corresponds to using 
one-half the value of the turbulent viscosity calculated from the TKE model 
in the k Jacobian, which was used for all three dependent variables. The 
matrix iteration was convergent, but the Jacobian distortion resulted in a 13% 
increase in CPU over reference case 1 and 24% increase in CPU over case 4. 
However, solution accuracy was not consequentially affected, as evidenced by 




table 7 

INFLUENCE OF THE JACOBIAN ON LINEAR FINITE ELEMENT SOLUTION “ BRADSHAW REUXING FLOW (TKE) 


Case 

No. 

Type of 
Jacobian 

CPU* 

Number 

of 

Passes 

Number of 
Iterations 
for 

First Pass 

AE/E 

AH/H 

1 

3 Jacobi ans 

1.00 

1291 

14 

.137 

-.059 

2 

u Jacobian 

1.08 

14S6 

14 

,151 

-.060 

3 

k Jacobian 

.91 

1236 

8 

.111 

-.058 

4 

. 0.5v* 

1.13 

1531 

12 

.116 

-.059 

5 

O.lv^ 

— 

- 

>30 

«M mm 



*NorniaT1 zed on case number 1 . . 
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comparing the normalized change in the energy norm for cases 3 and 4. Test 
case 5 corresponds to using only one-tenth the calculated turbulent viscosity 
within the Jacobian. This proved to be too inaccurate an evaluation, and 
convergence could not be achieved after iterating 30 times at the first 
integration step. Hence, a completely accurate evaluation of the Jacobian is 
not necessary to achieve an adequately-accurate engineering solution, and 
significant solution economies can result from taking advantage of the 
versatility embedded within the Newton iteration algorithm. 

Hyperbolic Equation Solution 
Equations Solved 

In Cartesian coordinates, the partial differential equation system 
governing transport of a scalar field, for example the transient continuity 
equation, is 

Uq) - ti'7q = 0 ■ (86) 

with boundary conditions 

A(q) = a^q + Vq • n + a^ =0 (37a) 

and an initial condition 

q(t.Q) - (87b) 

The goal of this analysis is a study of accuracy and solution economy of a 
factored Jacobian form of the developed Newton iteration-finite element 
solution algorithm. Select diver gen ce-fv^^ee rotational and irrotational 
velocity fields selected for this purpose include 

Constant: 

\ « Ujai + 0j) 


( 88 ) 
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Solid Body Rotation: 

iJg rfle 


Irrotational Flow About a Cylinder: 



Irrotational Cylinder Flow with Circulation: 





r^ 




(90) 


(91) 


In equations (88)-{91), is a reference freestream velocity, is the 
constant angular velocity, the two-dimensional solution spans 0 £x <^a, 

0 £y < b, is the circulation, and r is the cylinder radius. The initial 
distribution of q(x,y,0) is established as a "cosine-hill" rotated about its 
centre i dal node as 


qQ(x,y,0) = 100 



(92) 


where 0 £ r £ X is the local radial coordinate with origin (x^jy^), and X 
is spanned by M finite elements. Figure 44a illustrates the initial condition 
given by equation (92) for M=8 and (x ,y ) = (7,7) on a 32x32 uniform 
square mesh of span 0 £ a,b £ 80,000 m. 

The statement of the finite element solution algorithm for 
equations (86)-(87), on J2 = R^xt is, cf. Baker (1978), 


I EB200]{Q}g + [{U}^[B3001] + {V}g[B3Q023){Q}g = {0} (93) 


The Jacobian of the matrix Iterative solution of equation (93) is 




[Jl 


*-S|^[B20Q] + he|{U}^ [8300X3 +' {V}^[B30Q2l]j 
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(94) 


and the Iteration vector {SQ} is solved as 


[J]{5Q} = -{F} 


(95) 


For the factored Newton iteration algorithm, the dacobian [d] and 
elements of fF} are reexpressed on two-dimensional space in terms of the 
tensor matrix product (®). The two-dimensional factored Newton iteration 
algorithm is then written in the form 

(96) 

where Q represents an intermediate solution. In hypermatrix form, for a 
general one-step integration algorithm, equation (96) is written as 




Ai L4200] + h6A^ {YVCA30013 
®2 ®2 ® 


A^ [A2003 + h6A^ {UV [A30013 
®1 ®X ® 


mn 


pti 


d-t-l 


0 


Ag^rA2q03 - (Q>J * H[6yu>ItA300iKi5>P^j 

■J* Ci-e}A^ mJlASQGlKQjJ 

Ag MEOOl tOj+i] + ti. eAg {V}J,[A300n tQI? 


(1- e)Ag^{V}^[A30QlKq}^-+^^^^ 
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The comparisons between the factored Jacobian iteration algorithm, 
equation (96), and the conventional multi -dimensional algorithm, equations 
(93)-(94), are obtained for the velocity fields given in equations (88)-(89). 
Figure 44a shows the initial distribution of the wave packet on a 32x32 
uniform grid for velocity field u^^- Figure 44 b presents the conventional 
multi -dimensional bilinear finite element algorithm solution after 150 time 
steps with At =125 s. Figure 44c illustrates the final solution obtained by 
the multi -dimensional algorithm, but with the initial-value matrix diagonal- 
ized. Figure 44d shows the final solution- obtained with the factored Newton 
iteration algorithm, equation (96). It is virtually identical to the 
conventional results, and was obtained at approximately one-fourth expendi- 
ture of computer CPU and one-fifth the computer core requirement. These 
differences are essentially direct reflections of matrix bandwidth of [J] , 
hence become progressively more favorable as the mesh is refined. Table 8 
summarizes a comparison between the factored algorithm 'A', and the 
conventional multi -dimensional algorithm ‘B' for different values of the 
Courant number, uAt/A^, and two initial -value matrix structures. The CPU 
times for 'B' are five to seven times larger than for ‘A'. The finite 
element algorithm (S = 0) retains the peak better and has smaller trailing 
wakes than the diagonalized algorithm (S = 2). Numerical diffusion and 
dispersion error is also less for 'A' than for 'B' for largest value of 
Courant number. 

Corresponding accuracy, CPU and storage trends were obtained for 
velocity field hence, only the factored algorithm results are presented. 

-y. 

The solid-body rotation flowfield Ug is considerably more demanding, and 
provides* a quantization of dispersion error. The solution parameters (qg, a, 
b, At, M) remain identical, the diagonalized algorithm is relegated to 
history, and Figure 45 illustrates the solution obtained at the quarter. 


TABLE 8 


COMPARISON BETOEEN THE WCTORED NEWTON ITERATION ALGORITHM AND THE 
MULTI-DIMENSIOrWL BILINEAR ALGORITHM 





Peak Value {% of Original) 

Maximum Wake (55 of Original Peak) 

Courant 

Number 

CPU* 

A 

B 


A 

B 


. A 

B 

S » 0 S « 2 

S “ 0 

S » 2 

S « 0 S a 2 

S a 0 

S a 2 

0.1 

4.26 

19.52 

102. 66. 
(88.)^* 

102. 

62. 

(86.) 

3. -20. 

3. 

-19. 

0.5 

1.00 

4.5b 

102. 63. 

(87.) 

99, 

60. 

(84.) 

"5, -21. 

-5. 

-23. 

1.0 

0,52 

3.67 

90. 56. 

(92.) (82.) 

78. 

' (94,) 

53. 

(82.) 

-10. -21. . 

-19. 

-30, 


s § 


♦Normalized on FE factored solution for C«0.5, 
♦♦Maximum observed value. 
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three-quarter and full 360® rotation of the concentration packet. The 
initial -distribution would appear identical to Figure 44a, moved to the 
9 o'clock position, and the exact solution would be Lagrangian advection of 
the initial distribution without distortion. The factored algorithm is 
again essentially free of numerical diffusion (the peak level remains 
intact). The ripple structure in the ground plane is dispersion error, and 
while modest in comparison to other solution algoriths, cf. Long and Pepper 
(1976), is borderline on acceptability. Filters can be constructed to 
annihilate short period waves, e.g. Raymond and Garder (1976), and Figure 
44d illustrates the substantial imporvement accrued at the three-quarter 
turn of the filtered factored algorithm. Filtering can induce numerical 
diffusion, and the peak value has been reduced by 2%; a somewhat modified 
immediate trailing wake remains identifiable. For the multi -dimensional 
bilinear algorithm on a 32x32 uniform grid, the storage required for the 
Jacobian was 67xl089 locations compared to 2x3x1089 locations for the 
factored algorithm. CPU for the bilinear algorithm was ten times larger 
than for the factored algorithm, and there was essentially no difference 
between the two solution accuracies. 

The assessment of numerical preservation of symmetries and skew- 
symmetries by the factored algorithm was obtained using the irrotational 
velocity fields and Both correspond to flow about a cylinder of 
diameter 8 a^ and centered at the grid centroid. Two concentrations were 
symmetrically palced about the stagnation streamline, and all other solution 
parameters remain identical (q^^, a, b, U„, At, H). The computed filtered 
factored solution using 15^ is shown in Figure 46 for select time steps. The 
far- downstream peak values are within 2 % of the initial level, and dispersion 
error extrema are approximately ±4^. Exact symmetry of the two concentration 
cases throughout the advection was retained by the factored algorithm. 


56 


Figure 47 illustrates select results obtained by the factored algorithm, for 
the irrotational velocity field with circulation. Very large gradients 
are illustrated supported with acceptable dispersion error and negligible 
loss of peak value. 
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Figure 4. Slug Profile Used as Initial Condition 
for the Streamwise Velocity 
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Figure 7. Accuracy and Convergence with Discretization Refinement 
in Skin Friction Norm for a Parabolic Laminar Flow, 

Linear, Quadratic and Cubic Finite Elements 
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Accuracy and Convergence with Discretization Refinement 
in Energy Norm, laminar Boundary Layer, Linear Elements 
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Figure 12. Accuracy and Convergence with Discretization Refinement 
in Shape Factor Norm, Laminar Boundary Layer, 

Linear Elements 
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Figure 15. Accuracy and Convergence with Discretization Refinement 
(Nuiriber of Nodes) in Energy Norm, Laminar Boundary 
Layer, Linear, Quadratic and Cubic Finite Elements 
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Figure 16. 


Accuracy and Convergence with Discretization Refinement 
in Shape Factor Norm, Laminar Boundary Layer, Linear, 
Quadratic and Cubic Finite Elements 
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Figure 17. Accuracy and Convergence with Discretization Refinement 
In Skin Friction Norm, Laminar Boundary Layer, Linear, 
Quadratic and Cubic Finite Elements 
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Figure 20. 


Refinement 

in Skin Friction Norm, laminar Boundary Layer with 
ASeSly Linear Finite Elements, Consistent 
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Accuracy and Convergence with Discretization Refinement 
in Energy Norm, Laminar Boundary Layer with Pressure 
Gradient, Linear Elements, Diagonalized Initial-Value 





NormaMzed Error In Shape Factor Norm 


Ap X 10^ SI gn 

Go 

A 1.525 H^>H 

□ 3.0SS H^>H 

Non-Uniform 


A 


A< 


□ -y 


o . 


/ ■ / 

A ^ 

\ / ^ / 

N , \E + 

■ \: 
u' 


k=l 

S=0 


Discretization Refinement - tSJ 


Accuracy and Convergence with Discretization Refinement 
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Pressure Gradient, Linear Elements, Diagonalized 
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Figure 24. Accuracy and Convergence with Discretization Refinement 
in Snergy Norm, Laminar Boundary Layer with Pressure 
Gradient, Linear Elements, Finite Difference Initial- 
Value Matrix 
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Figure 28, Accuracy and Convergence with Discretization Refinement 
in Shape Factor Norm, Laminar Boundary Layer with 
Pressure Gradient, Quadratic Finite Elements 
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Figure 29. Accuracy and Convergence with Discretization Refinement 
in Skin Friction Norm, Laminar Boundary Layer with 
Pressure Gradients Quadratic Finite Elements 






Figure 30. Computed and Experimental Boundary Layer Parameters, 
Bradshaw Relaxing Flow, Linear Elements, MLT Closure 








Figure 32. Longitudinal Velocity Profiles, Bradshaw 
Relaxing Flow, Linear Finite Elements, 
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Figure 34. Finite Element Solution Energy* Bradshaw Relaxing Flow, 
Linear Elements, Consistent Assanbly, MLT Closure 
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Figure 36. Finite Element Solution Energy, Bradshaw Relaxing Flow, 
Linear and Quadratic Elements, Consistent Assonbly, MLT 
Closure 
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Figure 37. Accuracy and Convergence with Discretization Refinement 
in Energy Norm, Wieghardt Flat Plate Flow, Linear 
Elements, MLT Closure 
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Figure 40. Accuracy and Converg^^^^ Discretization Refinement 

in the Skin Fricti on Norm, ^ Plate Flow, 

Linear EiementSi 
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Accuracy and Convergence with Discretization Refinement 
in the Skin Friction Norm, Wieghardt Flat Plate Flow, 
Linear and Quadratic Elenients- MLT Closure 
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Figure 44. Advection of a Concentration Packet in 
Constant Velocity Field Ui, C = 0.1 









One-Quarter Turn 


nAt = 240 

Three-Quarter Turn 


nAt » 240 
Filtered 


Figure 45. Advectlon of a Concentration Packet in Solid 
Body Rotation Velocity Field U«, C = 0.1 
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Figure 46. Advection of a Concentration Packet Pair ^ oTiAiTra 
in Irrotational Velocity Field S 3 . C-O.l QWAWUB 
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Figure 47. Advectlon of a Concentration Packet Pair in Irrotational 
Velocity Field 0* with Circulation. C = 0.1 







